#include "NewtonInterpolation.h"

namespace NumericalAnalysis
{
    // Returns the f[x0,...,xk] used by divided differences.
	// Input: [x0, y0] ... [xk, yk]
    // Output[j] = f[x0, ..., xj]
	Utils::Matriz DividedDifference(const Utils::Matriz& pairs)
    {
    	Utils::Matriz result(pairs.mRows, 1);

		// result = [y0, ..., yk]
    	for(uint32_t i = 0; i < pairs.mRows; ++i)
    		result.mData[i] = pairs.mData[i * pairs.mColumns + 1];

    	for(uint32_t j = 1; j < pairs.mRows; ++j)
    	{
    		for(uint32_t i = j; i < pairs.mRows; ++i)
    		{
				// result[i] = f[xi+1, ..., xk] - f[xi, ..., xk-1] / xk - xi 
    			result.mData[i] = (result.mData[i * result.mColumns] - result.mData[(i - 1) * result.mColumns]) 
    				/ (pairs.mData[i * pairs.mColumns] - pairs.mData[(i - j) * pairs.mColumns]);
    		}
    	}

    	return result;
    }

    // Returns the newton productoria (x - xj) with 0 <= j < n
    double NewtonProductoria(const double x, const Utils::Matriz& pairs, const uint32_t n)
    {
    	double result = 1.0;

    	for(uint32_t i = 0; i <= n; ++i)
    		result *= (x - pairs.mData[i * pairs.mColumns]);

    	return result;
    }

    // Returns the pairs of (xi, Pn(xi)), where Pn 
    // is the Newton polynomial of interpolation.
    Utils::Matriz NewtonInterpolation(const std::string& pairsFilePath, const uint32_t n, 
                                      const std::string& valuesFilePath, const uint32_t m)
    {
    	const Utils::Matriz pairs = Utils::GetMatrixFromFile(pairsFilePath, n, 2);
    	const Utils::Matriz values = Utils::GetMatrixFromFile(valuesFilePath, m, 1);

    	Utils::Matriz results(m, 2);
    	double *resultsData = results.mData;

		// f[x0, ..., xn]
		Utils::Matriz dividedDifferences = DividedDifference(pairs);
    	for(uint32_t j = 0; j < values.mRows; ++j)
    	{
    		double currentResult = 0.0;
    		
			// Pn(xi) = f[x0] + f[x0, x1] (xi - x0) + ...
    		for(uint32_t i = 1; i < pairs.mRows; ++i) {

    			currentResult += dividedDifferences.mData[i * dividedDifferences.mColumns] 
    				* NewtonProductoria(values.mData[j * values.mColumns], pairs, i - 1);
			}

    		resultsData[0] = values.mData[j * values.mColumns]; // xi
    		resultsData[1] = currentResult; // Pn(xi)
    		resultsData += 2;
    	}

    	return results;
    }
}
